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Abstract 

Thermalization of quark-gluon plasmas in heavy-ion collisions is a difficult theoretical problem. 
One theoretical goal has been to understand the physics of thermalization in the relatively sim- 
plifying limit of arbitrarily high energy collisions, where the running coupling a s is weak. One of 
the current roadblocks to achieving this goal is lack of knowledge about the behavior of plasma 
instabilities when particle distributions are highly anisotropic. In particular, it has not been known 
how the magnetic fields generated by plasma instabilities scale with anisotropy. In this paper, we 
use numerical simulations in a first attempt to determine this scaling. 



I. INTRODUCTION AND RESULTS 



How do non-abelian plasmas that start far from equilibrium, such as quark-gluon plasmas 
produced in heavy ion collisions, equilibrate? This question has proven difficult to answer in 
detail even in the theoretical simplifying limit of weak coupling, appropriate to arbitrarily 
high energy collisions. A pathbreaking first attempt was made by Baier, Mueller, Schiff, and 
Son [1], who analyzed equilibration in the weak coupling limit via scattering processes of 
individual particles. The resulting picture of quark-gluon plasma equilibration is known as 
the bottom-up scenario. It was later realized, however, that collective effects in the form of 
magnetic plasma instabilities, known as Weibel or filamentary instabilities, necessarily play a 
role in bottom-up equilibration [2]. 1 These instabilities have long been known in traditional 
plasma physics [5], but their non-abelian counterparts develop somewhat differently. The 
effect of non-abelian interactions on the late-time development of plasma instabilities has 
been studied over the past few years with numerical simulations [6-13]. A great deal has been 
learned, but these simulations have not significantly explored the extremely non-equilibrium 
conditions relevant to the initial phase of bottom-up thermalization. 

In particular, Weibel instabilities are generated by anisotropic distributions of plasma 
particle momenta, as measured in local plasma rest frames. So far, simulations have mostly 
focused on the case of moderate anisotropy. 2 The bottom-up scenario, however, generates 
parametrically extreme anisotropies early on, before thermalization is achieved. Thermal- 
ization, of course, eventually produces isotropic (thermal) momentum distributions in local 
plasma rest frames. As an example, in the original bottom-up scenario (ignoring plasma 
instabilities), at one particular pre-thermalization moment of the expansion, the local dis- 
tribution of particle velocities looks like a pancake in momentum space, with 

Pz~9P±, (1-1) 

where g is the QCD coupling constant and z is the beam direction. Formally, in the limit 
of arbitrarily weak coupling g, this represents extreme anisotropy. To understand equilibra- 
tion in the weak coupling limit, one must therefore understand the development of plasma 
instabilities for the case of extreme anisotropy, p z /p± *C 1. The purpose of this paper is to 
make a first attempt to explore this limit using numerical simulations. 

Discussion of weak-coupling thermalization starts from the saturation picture of high- 
energy heavy ion collisions, where there is initially a non-perturbatively large phase-space 
density / ~ 1/g 2 of low x gluons with momentum of order the saturation scale Q s . These 
initial gluons are the "hard" particles in discussions of thermalization. We formally consider 
the case where Q s is so large that the running coupling a s (Q s ) can be treated as arbitrar- 
ily small. Bottom-up thermalization describes what happens as the plasma subsequently 
expands one-dimensionally between the two retreating pancakes of nuclear debris. The ex- 
pansion reduces the density of hard particles enough that one can treat them perturbatively 
for times r 1/Q S . In the first stage of the original bottom-up picture of thermalization, 



1 For a sample of earlier discussions of the possible role of Weibel instabilities in quark-gluon plasma 
thermalization, see Refs. [3, 4]. 

2 Two exceptions are the paper of Bodeker and Rummukaincn [11], with similar methods and aims to the 
current work, and the paper by Dumitru, Nara, and Strickland [10], which focuses on an initially perfectly 
planar distribution which is allowed to dynamically broaden with time. 
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which corresponded to 1 <C Q s r <C g~ 3 , the one-dimensional expansion effectively red-shifts 
the component p z of hard particle momenta along the beam axis, as measured in local 
plasma rest frames. For free particles, the expansion would drive the system away from 
local anisotropy as 

Pz 

— ~ (Qs'r) 1 - (free streaming) (1.2) 

V 

However, small- angle 2-^2 collisions between the hard particles tend to broaden p z /p, soft- 
ening the anisotropy to 

— ~ (<5 s t)~ 1//3 (original bottom-up) (1.3) 

in the original bottom- up analysis of Baier et al. [1]. This is a balance between one- 
dimensional expansion driving the system away from isotropy and collisions driving it to- 
wards isotropy. In Baier et a/.'s analysis, this relatively simple state of affairs continues until 
parametrically late times Q s r ~ g~ 3 , when other interesting things start to happen to bring 
about the eventual thermalization of the plasma. 

Plasma instabilities already play a role in the relatively simple first stage of bottom- 
up thermalization, however, and we will focus on this stage to motivate our investigation. 
In particular, plasma instabilities provide another mechanism to drive the system towards 
isotropy, and they change the exponent in (1.3). Weibel instabilities are associated with the 
creation of large, soft magnetic fields, which randomly bend the directions of the particles. 
How much bending occurs depends on the size of these magnetic fields B. Unfortunately, 
the parametric size of B in the case of extreme anisotropy (p z /p <C 1) has not been clear. As 
an example, there are two different guesses that have been made in the literature [14, 15], 
which would modify the original first-stage bottom- up behavior (1.3) to 

ft /(Qsr)- 1 / 4 , Ref. [14]; 

P l(Qsr)- 1/8 , Ref. [15]. 1 • ' 

The goal of this paper is to make a first attempt at resolving the issue by measuring the 
dependence of the soft magnetic fields B, caused by non-abelian Weibel instabilities, on the 
anisotropy of the hard particle distribution. 



A. Review: The limiting size of unstable magnetic fields 

Let fo(p) be the phase-space distribution of particles in the plasma, so that the density 
n is 

» = /(giA(p). (1.5) 

For moderately anisotropic fo(p), there is a single parametric scale of soft physics in the 
plasma which characterizes plasmon masses, Debye screening, and Weibel instabilities. For 
definiteness, we can take the scale of soft physics to be the effective mass of hard gluons 
in the plasma, given by [16, 17] 
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where there is an implicit sum over species, tR is a group factor, and v counts the number 
of non-color degrees of freedom (e.g. spin) for a given species. 3 For moderately anisotropic 
/o, the typical instability wavenumber /c U nstabic and growth rate 7 are both of order m^. 
Perturbation theory can be used to study the growth of instabilities from small seed fields. 
These instabilities cease to grow when their magnetic fields become large enough that their 
non-abelian self-interaction becomes important and perturbation theory breaks down [7, 
12]. Crudely speaking, that happens when gauge fields become important in soft covariant 
derivatives D = d — igA ~ i(k — gA), so that 



k 
9 



r\j — r^j 



(moderate anisotropy) (1.7) 



and 
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k 2 

B* ~ kA ~ — ~ 'jj^o (moderate anisotropy) (1.8) 

We write B with an asterisk subscript to denote, roughly speaking, the limiting size of the 
magnetic fields associated with unstable modes. This excludes other (higher momentum) 
modes which are excited at late times, associated with a cascade of plasmons that we will 
review later. 

For extremely anisotropic distributions fo(p), we have an additional parameter in the 
problem: the amount of anisotropy. Motivated by the application to bottom-up thermal- 
ization, we will focus on oblate distributions that are axi-symmetric about the beam axis z, 
and we will roughly characterize the amount of anisotropy by the typical magnitude 

0=^- = \v,\. (1.9) 
V 

Since 9 is parametrically small in the first stage of bottom-up thermalization, we need to 
know how the physics of instabilities depends parametrically on 6. A perturbative analysis 
of the instability shows that typical unstable modes have wave numbers k and growth rates 
7 of order [2] 

(k±,k z ) ~ (moo.tax) ~ ( m °°'~7r)' (1.10) 
7~moo, (1.11) 

where 

Aw ~ (1.12) 
is the maximum value of k for unstable modes. 4 



3 For a plasma of gluons, v = 2 and tn = 3. The v s in Rcf. [17] is this paper's v times the dimension of the 
particle's color representation. 

4 A way to remember this is as follows. The physical role of the scale moo in the context of instability growth 
is that 1 / moo is the time scale for currents to build up large enough to have important back-reaction on 
the fields. That 7 oc moo follows immediately Currents only build up if particles remain in a coherent 
region of single-sign field for this time scale. In time 1/moo, particles travel a transverse distance 1/moo, 
so k± ~ moo; but they only travel a z distance of ~ O/mr^, so k z ~ rrioo/9. This is illustrated in Fig. 1. 
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~ l/m„ 

FIG. 1: A schematic picture of a hard particle crossing a region of coherent magnetic field. 

What has been unclear is the size of the fields when unstable modes become non- 
perturbatively large and cease to grow. Here is a simple, hand- waving generalization of (1.8) 
which reproduces a conjecture made in Ref. [15]. Soft covariant derivatives D z = d z — igA z 
and D± = d± — igA± will become non-perturbative when 

k z k mSLX k± TfioQ . 
A z ~ — ~ and A± ~ — ~ , f 1.13 1 

9 9 9 9 

corresponding to magnetic fields 

B ± ~ (k ± A z or k z A ± ) ~ fcm ^ m °° (1.14) 



and 



2 

mi 



B z ~ k ± A ± ~ (1.15) 



The transverse fields dominate, with 

~ ~ — — . (l-lo) 

g Og 

Other arguments for the result (1.16) can be found in Ref. [15]. 5 In contrast, an earlier 
discussion by Ref. [14] assumed that B* ~ iTi^/g as in (1.8). One of our goals will be to 
distinguish these two possibilities using simulations to investigate the exponent v in 

(1-17) 

where 

{0, Ref. [14]; 
1, Ref. [15]; (1.18) 
2, Nielsen-Olesen limit. 

Here we show a third magnetic scale for comparison, the Nielsen-Olesen limit. It is associated 
with Nielsen-Olesen instabilities and is discussed in Appendix A. 

As in Refs. [14, 15], one can use (1.17) to parametrically determine how particle scattering 
from these fields broadens p z /p, and balance this against the one-dimensional expansion to 



5 See specifically Sec. V of Ref. [15]. Readers of other sections of Ref. [15] should be aware that some of the 
arguments there may be overly simplistic. In particular, see Ref. [18] for related discussion. 
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determine how p z /p scales with time. In a chaotic system, the coherence length of the 
unstable magnetic fields will be of order their wavelength, and so be of order 

11 19 

l± ~ - — ~ and l z ~ — ~ (1-19) 

in the transverse and z directions, as depicted qualitatively in Fig. 1. The particles have 
velocity (v±, v z ) ~ (1, 9) and will take time 5t ~ l/m^ to cross such a region. In that time, 
the magnetic force F will change the particle's p z by 5p z ~ F z t ~ gBJ± ~ gB^/rrioo. In 
time r, the particle will random walk through N ~ t/1± such changes, giving a total change 
of order 

A Pz ~ 5p 2 ~ (m^r) 1 / 2 ^ ~ . (1.20) 

Woo 

This will broaden the particle distribution to 

9 = * „ *Zl „ « r ) V2 . (1.21) 

p p o u Q s 



Solving self-consistently for 9, 



^ ( >^y' ( -> (122) 



Now we just need to know how moo depends on time. This was determined for the first 
stage of bottom-up thermalization by very simple arguments in the original work of Baier 
et al. [1] and remains unchanged in the presence of plasma instabilities. Comparing (1.5) 
and (1.6), one sees that m 2 ^ ~ g 2 n/p ~ g 2 n/Q s . Initially, at saturation, n ~ Q^/g 2 . In the 
first stage of bottom-up, there is no significant change in the number of hard particles, and 
so hard particle number density n dilutes from this initial value by the scale factor Q s r of 
one-dimensional expansion, so that n ~ Ql/ g 2 (Q s r). Putting everything together, 

moo ~ r-^Ql' 2 . (1.23) 

Inserting this into (1.22) produces the scaling (1.4) of p z /p with time quoted in the intro- 
duction. 

Table I A summarizes a variety of weak-coupling predictions for the first phase of the 
original bottom-up scenario [1] as well as its modification due to instabilities as conjectured 
in Ref. [15], corresponding to the limiting field (1.16) above. Here, we have defined the 
dimensionless time f = Q s t. Since many readers may be more familiar with equilibrium 
plasma physics than with the scales of bottom-up thermalization, we also show, for the 
sake of qualitative comparison, what similar scales would be for (i) an equilibrium plasma at 
temperature T, and (ii) a "squashed" equilibrium plasma which has the same density n ~ T 3 
and typical energy p ~ T of particles but has particle momenta distributed anisotropically 
with p z /p I. In the thermal case, the hierarchy of different mass scales is controlled 
by the small parameter g. In the bottom-up case, it is instead controlled by the small 
parameter (Q s t) _1 - Note that the bottom- up scales satisfy the hierarchy that the instability 
growth rate is parametrically faster than both the expansion rate 6 and the rate for individual 



6 See Ref. [19] for a recent analysis of the unfavorable effects of expansion on instabilities for heavy ion 
collisions at realistic (rather than arbitrarily large) energies. 
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general 


thermal squashed 
thermal 


original guess 
bottom-up Ref. [15] 


hard particle momenta p 
particle isotropy 9=p z /p 
hard particle density n 
phase space density / 
hard plasmon mass TOqo 
particle collision rate 
expansion rate 


n/9p 3 

/ 2^ /„ 
Vfi ' n /P 

5 4 n(l+/)/m^ 


T T 
1 0< 1 

1 0- 1 


Q s no change 
f-V3 f-1/8 

r Qs/g no change 

f-2/3/^2 f-7/8/^2 

r ' Vs no change 
f-^Q s f- 7 /8Q s 
r _1 no change 


instability wave number 
instability growth rate 


moo W 




f- 1 /6Q s f -3/8Q g 
f _1 / 2 Q S no change 



TABLE I: A table of the parametric dependence of various scales for (i) a thermal distribution, (ii) 
a "squashed" thermal distribution with the same density n but extreme momentum anisotropy 9, 
(hi) the first stage (1 <C Q s t <C g~ 3 ) of the original bottom-up thermalization scenario of Baier et 
al. [1], and (iv) the changes to bottom-up due to instabilities based on the conjectured dynamics 
of Ref. [15]. In this table, f = Q s t and the phase space density / refers to the largest values 
of f(p) (and not to the angular- averaged values). The "particle collision rate" refers to the rate 
of individual (incoherent), small-angle 2^2 scattering of hard particles from each other. The 
instability wave number refers to k ~ k z . In contrast, k± ~ moo as discussed in the text. An 
asterisk (*) indicates general formulas that apply only to moderate to extreme anisotropy and not 
to isotropic or nearly isotropic situations. 



(incoherent) 2—^2 hard particle collisions. For the purpose of determining the limiting of 
(1.17) in the specific context of bottom-up thermalization in the weak coupling limit, this 
hierarchy allows one to ignore the effects of both expansion and individual collisions when 
simulating plasma instabilities [2]. 



B. Overview of simulation method and what we measure 

In order to cleanly separate scales in the weak-coupling limit, simulations are carried 
out for the hard-loop effective theory of soft excitations as in Refs. [6, 7, 12, 20]. This 
effective theory is a non-abelian version of the linearized Vlasov equations of traditional 
plasma physics, which are based on collisionless kinetic theory for hard particles coupled to 
a soft, classical gauge field. We use the formulation of Ref. [13], where the equations are 

D v F» v {x,t) = I v»W(v,x,t) , (1.24a) 

J V 

(D t + v D X )W = ml [E ■ (2v -V v )+B-(vx V v )] Q(v) . (1.24b) 

Here, the field W a (v, x) represents the net (adjoint) color of all particles moving in direction 
v at point x. The first equation is the Yang-Mills field equation, with W(v) giving rise 
to a current. The second equation, derived in this form in Ref. [12], shows how electric 
and magnetic fields can polarize the colorless distribution of particles to create a net color 
moving in each direction. In this equation, Q(v) is a static quantity which parametrizes the 
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angular distribution of the initial, background distribution fo(p) of hard particles. W(v, x, t) 
is generated by small fluctuations of f(p) from fo(p). (For the particular weak-coupling 
questions treated here, it is allowable to treat the hard particle fluctuations as small.) The 
dynamics of the soft fields is equivalent to that of hard-loop effective theory [20] . The use 
of classical equations (1.24) can be justified for the applications at hand because deBroglie 
wavelengths of the hard particles are parametrically small compared to the soft physics 
distance scales, and because the instability causes the soft gauge fields to grow parametrically 
large enough to be classical. 

To discretize these equations for simulation, we follow the methods of Refs. [12, 13], with 
a small but important change discussed in Sec. Ill to allow us to more efficiently simulate the 
case of extremely anisotropic distributions. Also, like previous studies of Weibel instabilities, 
all of our simulations will be for SU(2) gauge theory for reasons of computational simplicity. 
We expect this to be qualitatively similar to SU(3) gauge theory; we are not aware of any 
reason why they would be different. 

Fig. 2 shows an example from Ref. [12], showing the total energy density in soft magnetic 
fields as a function of time. This particular simulation was for moderate anisotropy and 
started from tiny initial conditions for the gauge fields. There is exponential growth at early 
times, due to the instability, and linear growth at late times. The linear growth does not 
represent continued growth of the unstable modes. Instead, the unstable modes stop growing 
but, through interactions, pump energy into a cascade of increasingly higher momentum, 
stable modes [13]. This cascade takes the form of a gas of plasma excitations of the classical 
gauge field with momentum q > ^unstable, which are perturbative for q ^> A: uns tabic- At late 
times, the total classical magnetic field energy density £^ ot = ^B 2 is dominated by the 
energy of these perturbative plasma excitations, rather than the energy density Ef ~ \B 2 
of the softer ik ~ /^unstable) unstable modes. For this reason, we cannot simply measure the 
total magnetic energy density |£> 2 at late times and take a square root to find the limiting 
size £?* of the magnetic fields associated with unstable modes. And it is the soft fields 
k ~ /cunstabie, not the higher momentum plasmon excitations, which dominate the scattering 
of hard particles and so determine the evolution (1.17) relevant to bottom-up thermalization 
[15]: Even though the k ~ instable modes carry less energy, they are more effective at 
scattering. To determine what we really want to know, we need some measurement other 
than the total magnetic energy density at a single late time. 

In this paper, we will use an indirect method to investigate anisotropy dependence which 
is relatively easy to implement. We will measure the slope dS^/dt of the late-time linear 
growth of the total magnetic energy density and determine how it scales with anisotropy. 
Here is a model of how one might expect this slope to behave. The source of increasing total 
magnetic energy comes from the unstable modes, which take energy from the hard particles 
and, through interactions, dump it into the cascade of plasmons. As a thought experiment, 
imagine that half the energy density £f in the unstable modes were abruptly transferred to 
the cascade of plasmons. How long would it take the unstable modes to grow back to their 
original, limiting size? Parametrically, the time should be of order the inverse instability 
growth rate t ~ 7 _1 ~ mZo- (Even though this is a perturbative estimate of the growth 
rate, it should still be parametrically correct in the region where perturbation theory starts 
to break down.) So the rate energy is pumped into the cascade can be expected to be of 



order 7£f : 



d£ B 



'tot _ 



d 

dt 



(^ 2 )»~^. B ~7B, 2 
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5 



oo 




dt 
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FIG. 2: Magnetic energy vs. time for a sample simulation from Ref. [12] for moderate anisotropy, 
starting from a very small seed for the Weibel instability. The figures are the same except that the 
vertical axis is logarithmic in the left-hand figure and linear in the right-hand figure. The solid line 
is the result for non-abelian gauge theory in three spatial dimensions. For comparison, the dashed 
line shows a simulation in an abelian theory, and the dotted line shows a non-abelian simulation 
in one spatial dimension. 

where we have used the parametrization (1.17) of 5*. By measuring how this slope scales 
with 9, we can extract the desired exponent v that determines £>*, assuming that the physical 
argument for (1.25) is correct. 

Following Ref. [13], we will generally start our simulations with large initial gauge fields 
so that we can quickly and easily get to the late-time limiting behavior. For extremely 
anisotropic hard particle distributions, this is a non-trivial choice: recent simulations [11] 
starting instead from tiny initial gauge fields find qualitatively different behavior. We will 
return to this point in Sec. II, where we argue that large initial conditions are appropriate 
to understanding how bottom-up thermalization is modified by instabilities. 

We should note that the limiting field -B* we have used to present a qualitative picture of 
the physics of instabilities is not a precisely defined quantity. Unlike the total magnetic field 
-Stotj we know of no unique, convention-independent, gauge-invariant definition of the mag- 
nitude of -B*, and so is only useful for parametric estimates. In contrast, the observable 
(1.25) discussed above is gauge- invariant. 



C. Overview of Results 

Fig. 3 shows an example of the total magnetic energy density £^ ot vs. time for three 
different anisotropies, starting from strong, non-perturbative initial conditions which we 
shall detail later. One way to parametrize the amount of anisotropy is to rewrite (1.12) as 

(1.26) 

"'max 

where we compute moo/A; max perturbatively for each background hard particle velocity dis- 
tribution VL{v) [i.e. each distribution fo(p)] that we simulate. Fig. 3 shows increasing slope 
for increasing anisotropy. We will then rewrite the scaling form (1.25) in the form 

g 2 dSg t /dt 

- constant, (1.27) 



" l oo ^max I* 
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FIG. 3: Examples of the linear growth of total magnetic energy with time. From bottom to top, 
the three curves have k raa:x _/m 00 = 2.16, 3.15, and 4.10 respectively (corresponding to the Nq = 3, 
5, and 7 distributions described in Sec. Ill A). 



where we take 7* to be the largest unstable mode growth rate computed in perturbation 
theory. (This rate approaches rrioo/^/2 in the limit of extreme anisotropy [2], but we have 
chosen to keep 7* explicit in our formula because the approach to this limit is a bit slow. 
Details will be given in Sec. Ill A.) 

Fig. 4 shows the left-hand side of (1.27) vs. our measure m^/km^ of anisotropy for a 
variety of different simulations, plotted with exponents v = |, 1, and §. Each point has 
systematic errors of order 15%. The v — 1 version plausibly approaches a constant in the 
extreme anisotropy limit m^/k^^ — > 0. The v = | and v = | figures clearly rule out v < ^ 
and v > |. Of the three possibilities v = 0, 1, and 2 considered in (1.18), we conclude that 
only v = 1 is consistent with this measurement. 



II. LARGE VS. SMALL INITIAL CONDITIONS 



We have initialized our simulations with large initial gauge fields. In contrast, many 
simulations in the past, such as Fig. 2, have started from tiny initial gauge fields in or- 
der to observe the crossover from perturbative, exponential growth of instabilities to the 
limiting late-time behavior. Fig. 5 shows the difference for moderate anisotropy: we have 
superposed the tiny initial condition simulation of Fig. 2 (Ref. [12]) with an otherwise iden- 
tical simulation starting from large intitial conditions. For tiny initial conditions, there is a 
significant spurt of continued exponential-like growth even after the field strength reaches 
non-perturbatively large values. It is only later, at much higher energy, that linear growth 
finally sets in. Bodeker and Rummukainen [11] have found that this spurt of post-non- 
perturbative exponential growth for tiny initial conditions becomes much more significant for 
extreme anisotropy. In their simulations for extreme anisotropy, they see only exponential- 
like growth at late times; they do not see late-time linear behavior at all. It is possible 
that the late-time behavior is ultimately linear but sets in at such large field energy that 
their simulations cannot reproduce it because of lattice spacing artifacts. But regardless, 
the full story of the development of plasma instabilities appears to be qualitatively different 
depending on whether or not one starts with large or tiny initial conditions. 

Which type of initial condition is relevant to a new scenario of bottom-up thermalization? 
We argue that it is the case of non-perturbatively large initial conditions that we have inves- 
tigated in this paper. Consider some time T\ in the first stage of a new bottom-up scenario 
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FIG. 4: The slope d£^ ot /dt of the linear growth of total magnetic field energy, measured in units of 
^maxT*/^ 2 ; as a function of anisotropy for (a) v = ^, (b) v = 1, and (c) v = §. Simulation 



4-2uh,2u 



parameters are listed in Table III in Section III; squares are the default values and crosses are 
alternate values at the bottom of the table. The numbers by the data points indicate the order 
Nq of the distribution, as described in Sec. Ill A. 



that accounts for instabilities, with Q s t\ ^> 1. We've already reviewed how the instability 
growth rate is large compared to the expansion rate 1/ti, and so the unstable modes will 
have grown to become non-perturbatively big (or perhaps larger). This population of unsta- 
ble modes is depicted very crudely by the curve in the cartoon of Fig. 6. (In addition, higher 
momentum modes may be populated due to having been unstable at earlier times, or due to 
interactions.) It's important to note that the unstable modes will grow to non-perturbative 
size regardless of how small the initial seed fields for those unstable modes are, because the 
instability growth rate is fast and quantum fluctuations will seed the unstable modes even 
if nothing else does. 7 

Now consider what happens a little later, at time r 2 = 2t\. The set of unstable modes 



7 In more detail, quantum fluctuations by themselves would correspond to fluctuations of order A ~ k max in 
the size of typical unstable modes. The instability would cause these to grow to non-perturbative size A ~ 
fcmax/s in time of order 7 -1 times the log of the size ratio: 7 _1 ln(l/g) <~ m^} ln(l/p) ~ (t/Q s ) 1//2 ln(l/</). 
This time is much shorter than the life r of the system when Q s t ^> ln 2 (l/(/), which we can roughly think 
of as the condition Q s t ^> 1 for the applicability of bottom-up thcrmalization, since we have generally 
not tried to keep track of logarithms in discussions of scales. 
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FIG. 5: The difference between tiny initial conditions, as in Fig. 2, and large initial conditions, 
similar to Ref. [13]. (The large initial conditions were set as in Sec. IIIC but with T = Am^/Sg 2 , 
^smear = ^fnooi and squeeze s = 1.5.) 



/ 




k 



t 

FIG. 6: A cartoon of the occupation number f(k) of unstable modes at time r±, with the shaded 
area depicting those modes which are perturbatively unstable at the later time r 2 = 2t\ . 

shrinks a bit in k space. Specifically, combining (1.12), (1.22), and (1.23), we have 

fc m ax~^~gs(g s r)-( 1+2 ^ 1+ ^. (2.1) 

For definiteness, consider v — 1, for which fc max ~ Q s (Qst) -3 ^ 8 - Then k max decreases by 
a factor of 2~ 3 / 8 ~ 0.77 when time r increases by a factor of 2, The unstable modes at 
the later time T2 are therefore shown by the shaded area in Fig. 6, and we see that they 
are already initialized with non-perturbatively large fields because of the earlier instability 
growth at time T\. 

We have tried to make our argument general: Regardless of whether we start with large 
or small seeds for instability growth at r = Ti, we will then get large seeds for instability 
growth at later times such as r 2 = 2t\. But the same argument means that we had large 
seeds at T\ also because of yet earlier instability growth at time 7i/2. We can follow this 
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argument all the way back to times of order the saturation time, when all the relevant modes 
of the fields started non-perturbatively large. We conclude that the typical unstable modes 
have high occupancies at all times during the initial stage of the new bottom-up scenario. 



III. SIMULATING EXTREME ANISOTROPY 

To simulate the non-abelian Vlasov equations (1.24), we need to discretize the arguments 
x and v of the fields W(x, v, t) and A(x, t). For x, we put the system on a spatial lattice. For 
velocity v, we follow Refs. [12, 13, 21] and expand in spherical harmonics Yi m (v), truncating 
the expansion at some maximum value £ max of £: 



L rnax 



W(x, v, t) = ^ ^ W tm {x, t) Y em (v), (3.1) 

£=0 m 

where our convention is to normalize the spherical harmonics so that the angular average of 
Y^ m (v) Y£' m '(v) is 8w8 mrn i. (We place the hat over Y^ m as a reminder of this non-standard 
normalization convention.) The axi-symmetric, hard particle background velocity distribu- 
tion Q(v) has the expansion 

^max ^raax 

n(v) = J2^eY m (v) = J2^£+l) 1 / 2 P e (v z ), (3.2) 

£=0 £=0 

where the Pe(x) are Legendre polynomials. The explicit form of the equations of motion 
(1.24) in terms of the W^ m 's is given in Ref. [12]. 

In practice, we must choose £ max large enough to obtain results close to the £ max — > oo 
limit. More anisotropic distributions Q(v) will require larger £ max and therefore greater 
computational resources (both memory and time, to store and evolve more W^ m 's). Below, 
we first describe our choice of distributions Q(v) to simulate. Then we explain and justify 
our method for making simulations of very anisotropic distributions practical, which is to 
reduce the number of W^ m 's by limiting m to \m\ < m max with m max ~ 6. 



A. Choice of hard particle distribution fo(p) 

For a given maximum £, we would like to find a velocity distribution fl(v) of hard particles 
which is as anisotropic {9 ~ \v z \ 1) as possible. To be physical, Q(v) should be non- 
negative. 8 After some experimentation, we settled on the following form, parametrized by 
an integer order N^: 

Of. ^ = - V ^ rE=i(«i - ^\ for N n = 2n+ 1; 

1 Z) \W ~ vlf nr=i(«, - for N Q = 2n + 2; 1 ' } 

where the normalization M is chosen to satisfy our convention that the angular average Qq 
of Q(v) is one. We choose the at to minimize (vf), performing the minimization numerically 
for each Nq. 



We do not know if there would be any problem for simulations if Q(v z ) had tiny negative values for some 
v~, but it seems safer to avoid this. 
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1 i 


\^ z j rms 








1 


0.4472 


l 


0.500 


0.111 


2 


0.3780 


1.414 


0.649 


0.191 


3 


0.2852 


2.155 


0.875 


0.310 


4 


0.2506 


2.542 


0.979 


0.354 


5 


0.2093 


3.149 


1.130 


0.408 


6 


0.1887 


3.542 


1.221 


0.435 


7 


0.1653 


4.099 


1.342 


0.466 


8 


0.1516 


4.500 


1.424 


0.484 


9 


0.1366 


5.031 


1.527 


0.505 


10 


0.1269 


5.438 


1.602 


0.518 


11 


0.1163 


5.954 


1.693 


0.533 


12 


0.1091 


6.365 


1.763 


0.543 


13 


0.1013 


6.870 


1.846 


0.554 


14 


0.0957 


7.285 


1.911 


0.562 


15 


0.0897 


7.783 


1.987 


0.570 


17 


0.0805 


8.693 


2.119 


0.583 


19 


0.0731 


9.601 


2.245 


0.594 


21 


0.0668 


10.508 


2.361 


0.603 


23 


0.0616 


11.414 


2.475 


0.611 


25 


0.0571 


12.319 


2.583 


0.618 



TABLE II: For each of the hard particle distributions designated by Nq, the quantity (v z ) rms = 
(v 2 ) 1 ^ 2 measures the narrowness of the velocity distribution about the transverse plane. k max is 
the maximum unstable wavenumber. 7* is the largest perturbative growth rate of the field modes 
A(k) and corresponds to wavenumber fc*. The corresponding perturbative growth rate of magnetic 
energy is 27*. 

The expansion of Q(v) involves spherical harmonics with I < Lq = 2Nq, and the cor- 
responding coefficients Qi are listed in Appendix B for various choices of Nq. We will see 
later, looking at the £ max dependence of results, that we can get close to the £ max — > 00 limit 
for a given N n using £ max > 1.5 L n = 3N n . 

Table II summarizes basic properties of these distributions for various values of Nq. 
Increasing anisotropy is signaled by decreasing (f 2 ) rms = (v^) 1 ^ 2 and increasing /c m ax/^oo- 
For graphical comparison, Fig. 7 shows, for different values of Nn, the perturbative growth 
rates of unstable modes as a function of wavenumber k in the case that k points exactly along 
the beam direction. For each distribution, /c max denotes the largest unstable momentum, 
7* is the largest growth rate, and is the corresponding momentum. Figs. 8 and 9 show 
kmax, lA^rmsj and 7* vs. Nq. For comparison, the moderately anisotropic distribution 
previously simulated in Refs. [12, 13] is roughly comparable to our Nq = 3 distribution, and 
the most extremely anisotropic distribution simulated in Ref. [11] is roughly comparable to 
our Nq = 15 distribution. 9 



Specifically, the distribution used in Refs. [12, 13] has (v z ) lms = 0.312, and the L asym = 28 distribution 
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FIG. 7: Perturbative instability growth rates ^(k) vs. k/rrioo for the values of Nq listed in Table 
II. Solid (dashed) black lines are odd (even) Nq, staring with Nq = 1 at the bottom and running 
up to Nq = 25 at the top. The horizontal dotted line is the maximum possible 7, which is 1/V2, 
and the dashed line approaching it is the case Nq = 00, given in Ref. [2]. 
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FIG. 8: The values of l/(f 2 ) rm s (circles) and k^^jra^ (crosses) plotted vs. Nq. 
B. A reduced set of Y^m's 

The expansion of extremely anisotropic distributions fl(v z ) in spherical harmonics Yi m (v) 
requires large i values but, due to the axial symmetry of the distribution, only the m value 
m = 0. The dynamics (1.24b) of the fluctuations W(v, x, t) in the distribution, however, will 
create M^ m 's with non-zero values of m. We might hope that only small m values turn out 
to be significant. Though there will be a lot of rapid variation in how W(v, x, t) depends on 
v z , because we are studying the case of extreme anisotropy, there might be relatively smooth 
dependence on (v x ,v y ). We will verify this picture below using simulation data, and also 
give some qualitative arguments why one might expect it. We can take advantage of this 
smooth dependence by placing an upper bound \m\ < m max on the range of m we include 



of Ref. [11] has (v z ) Tms = r]/\/S = 0.0864. Also, the notation m^/m§ of Ref. [11] is equivalent to our 
notation Qi. 
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FIG. 9: The maximum growth rate 7* (in units of m^) plotted vs. Nq. The top of the graph 
represents the Nq — > 00 limit of 7* — > moo/ a/2 [2]. 

in our simulations, so that the expansion (3.1) of W becomes 

^max 

W(x,v,t) = J2 W tm (x,t)Y em (v), (3.4) 

e=o \ m \<e 

\m\<m max 

For unrestricted m's, the total number of W tm at each lattice site (and so the resources 
required for the simulations) would grow quadratically with £ max . For a fixed bound 
\m\ < m max , however, they only grow linearly for large £ max , making simulations of ex- 
treme anisotropy practical. 

Fig. 10 shows an example of linear growth of total magnetic field energy for Nq = 7 
simulations with £ max = 24 and several different values of m max . As can be seen, m max = 6 
is large enough to reproduce the correct (m max — > 00) slope, and this is the value of m max 
we will use in our simulations. 

In previous work on simulations for moderate anisotropy [13], we found that the system- 
atic errors arising from a finite cut off £ max on £ could be reduced by damping the dynamics 
of modes with I near the cut-off. We have slightly improved this method and extended it to 
apply also to the new cut-off m max on m. Details are given in Appendix C. Such damping 
has been used in all the simulations reported in this paper. 

The real test of the viability of using relatively small m max cut-offs comes from simulations, 
such as Fig. 10. However, one can get some rough idea of why it can work by considering 
perturbative formulas for some of the important features of unstable modes and the resulting 
cascade of plasmons. If we treat the gauge field perturbatively in the W equation (1.24b), 
replacing by d^, we can Fourier transform from (x,t) to (k,u) and then solve for W: 

W(v, k, uj) = im 2 ^ (u-v- k)' 1 [E(k, u) ■ (2v - V„) + B(k, u) ■ {v x V„)] n(v g ) . (3.5) 

Together, the factors to the right of the (u — v ■ k)^ 1 in this formula only generate v± 
dependence with |m| < 1. All higher m components in the result for W are generated by 
the factor 

(uo — v ■ k)^ 1 = (uo — Vj_ ■ k± — Vzkz) -1 . (3.6) 

Now consider the dominant unstable mode. As mentioned earlier, this mode has k along 
the z axis (for the type of anisotropy we consider), and so k± = 0. Then the factor (3.6) 
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FIG. 10: Linear growth of the total magnetic energy with time for several different values of r^max • 
The parameters are the same as our canonical Nq = 7 simulation except the box size is 40 2 x 32 
and the squeeze factor is only 2. 



has no v± dependence, and so the W field which describes the dominant instability involves 
only \m\ < 1. 

As another example, consider the dispersion relation of transverse plasmons. A standard 
method for deriving the dispersion relation is to insert the result for W into the Yang-Mills 
equation (1.24a), which generates the hard-loop self-energy correction to the vacuum relation 
uj 2 = k 2 . How much will we disturb this dispersion relation if we throw away modes of W 
with m > m max ? For high momentum plasmons (k ^> moo), such as those that dominate the 
cascade at late times, the effect is tiny simply because all medium effects to the dispersion 
relation are tiny in this limit. For very low momentum plasmons (k <C uj ~ moo), we can 
ignore the v ■ k altogether in (3.6), and then the W field will again have only |m| < 1 
components. It is only for intermediate momentum plasmons (k ~ moo) that finite m max 
does violence to the plasmon dispersion relation. However, (3.6) is a fairly smooth function 
of v± in this regime because the denominator never gets close to zero for k ~ moo plasmons 
(u; and to — k are both of order moo), and so (3.6) and therefore W can be reasonably 
approximated by a superposition of relatively low m's. 

Finally, we should check that we have chosen large enough values of £ m3jX in our simula- 
tions. In general, we find that £ max ~ 3Nq is quite adequate. As an example, Fig. 11 shows 
the £ max dependence of Nq = 7 simulations for fixed m max = 6. Our standard simulation 
choice for £ max is 24 for Nq = 7. See Table III for our default simulation parameters in other 
cases. 



C. Initial conditions 

Following Ref. [13], we use strong, non-perturbative initial conditions for the magnetic 
field B, so that the system starts linear energy growth behavior as quickly as possible. The 
electric and W fields are, for simplicity, initialized to zero. In order to see the linear energy 
growth associated with cascade development as early as possible, it is advantageous to choose 
initial conditions which do not significantly populate modes with large wavenumber. 
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FIG. 11: Linear growth of the total magnetic energy with time for several different values of £] 
and fixed m max = 6, for the hard particle distribution Nq = 7. 
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32 
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64 2 x 32 


11 
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3.5 


64 2 x 32 


13 
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64 2 x 32 


15 


1.0 


56 
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25 
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80 
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64 2 x 28 
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24 


2.5 


64 2 x 32 
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24 
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24 
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64 2 x 32 
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24 
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64 2 x 32 



TABLE III: The default parameters for our simulations and their initialization, as a function of 
Nq. The corresponding values of fc max are given in Table II. Other default parameters include 
m max = 6, initial temperature T = k max /g 2 , and initial smearing wavenumber k srneajT = k max . The 
simulations below the horizontal line correspond to the crosses in Fig. 4. 

In the moderate anisotropy simulations of Ref. [13], the initial magnetic field was con- 
structed by taking a thermal initial state with temperature T = 2m 00 /g 2 and then perform- 
ing gauge-invariant smearing (sometimes called cooling) of the configuration to eliminate 
wavenumbers k 3> m^. In perturbative language, this cooling corresponds to replacing the 
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initial thermal field A t h e rm by 



A(k) = A thcrm (k) exp(-/c 2 /A; s 2 mear ), (3-7) 

where r = ^/k 2 mear is the smearing parameter. In Ref. [13], we chose k smeair = Im^. 

Here, we follow a similar procedure, but the unstable modes that we want to initially 
populate are generally more extremely anisotropic, having (k±, k z ) ~ (moo, fc ma x) with 
fc ma x 3> ^oo- We have found that it helps to arrange a related anisotropy of our initial 
fields by squeezing the initial distribution in the z direction. In perturbative language, our 
initial choice corresponds to 

(A ± , A z ) [k] = ( A ±jthcrm / s, A zMlcrm ) [sk±, k z ] exp ( - (s 2 k]_ + k 2 z )/k 2 mc ^ , (3.8) 

where s is the squeezing factor. 10 In our simulations, we have generally chosen T = k max /g 2 , 
^smear = k max , and s between 1.5 and 3.5 depending on the amount of anisotropy. See Table 
III for our default simulation parameters. 

D. Lattice spacing and volume 

It is important to check that the lattice volume is large enough to be in the infinite volume 
limit and the spacing is small enough to be in the continuum limit; otherwise the lattice 
calculation is not simulating the desired continuum physics. It would be prohibitive to check 
this at every lattice spacing, so we have "spot checked" this at a few levels of anisotropy, 
with the most thorough study at N n = 7 and N n = 15. 

For highly anisotropic lattices, the physical scales possibly relevant to out problem para- 
metrically span a range from to k max . One might worry that, as particle distributions 
are taken more and more anisotropic, it becomes harder and harder to span these scales 
with a computationally practical lattice. Naively, to be perfectly safe, we would like physi- 
cal lattice dimensions L ^> 27r/m 00 and lattice spacings a <C 2n/k max . In this section, we'll 
see how well we do with lattices of practical size. 

1. Physical volume 

Figs. 12 and 13 show the volume dependence, at fixed lattice spacing, of the evolution of 
magnetic energy with time. The first figure is for the N n =7 hard particle distribution. The 
second figure is for Nq—15, the second most anisotropic distribution included in our results 
of Fig. 4. Our default lattice size of 64 2 x 32 corresponds to approximately (15.6/moo) 2 x 
(7.8/raoo) for Nq = 7 and (8.2/moo) 2 x (4.1/toqo) for Nq = 15. For small volumes, the 
simulations produced exponential rather than linear growth. 11 But linear growth appears at 



Our technical procedure is to choose the initial magnetic field by the usual procedure but pretending that 
the lattice is asymmetric with lattice spacing a± = a/s in the transverse directions, compared to a along 
the z axis. We then re- interpret the resulting initial condition as living on the symmetric lattice (a z = a±) 
used in our simulations. 

These are small volumes with periodic boundary conditions. One should not expect this small- volume ex- 
ponential growth behavior for a comparably small volume of hard particles in infinite space, surrounded by 
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N n =7 L„,m ma =24,6 squeeze = 2(3 for 64 2 ) k m „ v a=l 

u max' max ' i \ / max 




FIG. 12: Linear growth of the total magnetic energy with time for several different physical volumes, 
at fixed lattice spacing, for the hard particle distribution iVr> = 7. 
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FIG. 13: As Fig. 12 but for the more anisotropic distribution Nq = 15. 



large enough volume, and a comparison of the large volume curves suggests that our default 
lattice size of 64 2 x 32 is adequate, even for our highly anisotropic distributions. 

In order to be able to run our simulations on desktop computers, we by default took the 
physical lattice size L z in the z direction to be half that in the x and y directions. This 
choice is motivated by the fact that, in the highly anisotropic case, unstable modes have 
parametrically smaller wavelength in the z direction (~ l//c max ) than in the perpendicular 
directions (~ 1/rrioo). Of course, that doesn't exclude the possibility that stable modes with 



vacuum. In that case the hard particles would escape the small volume within the time scale characteristic 
of the instability growth. 
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FIG. 14: Linear growth of the total magnetic energy with time on an L\ x L z lattice for several 
different choices of L z (at fixed lattice spacing), for the hard particle distribution Nq = 7. 



size k x ~ k y ~ k z ~ l/m^ might be important in the development of linear growth, and so 
we should investigate the matter with simulations. Fig. 14 isolates the effect of varying L z 
while holding the L x and L y fixed. (Note that L x = L y is smaller here than in Fig. 12.) An 
L z that is half of L x = L y appears adequate for reproducing the large- linear slope. 



2. Lattice spacing 

Fig. 15 shows how our simulations depend on lattice spacing for fixed physical volume, 
for the distributions Nq = 7 and A*q = 15. In order to isolate the effect of lattice spacing, we 
have used the same initial conditions for all these simulations. More precisely, we generated 
initial conditions for the finest lattice (96x96x48), and then we used blocking to generate 
similar initial conditions for the coarser lattices. 12 

At all but the finest lattice spacing in Fig. 15a, one can see some curvature to the late-time 
"linear" growth behavior. This curvature is a lattice artifact, but it means that we need a 
procedure for extracting a single "slope" from such simulations, in order to present results 
such as Fig. 4. Note that the coarsest lattice spacing results in Fig. 15 look like sections of 
tanh curves, after an initial transient. Inspired by this observation, we have chosen to fit 
each of our energy curves to the form 

\B 2 (t) = atitanh ( ^— ^ J (3.9) 



For instance, to block by a factor of 2 in every direction, one could replace appropriate pairs U\ and U% of 
consecutive links by a single link U\Ui. In practice, we use the slightly improved method of averaging this 
with the four "staples" that move one link transversely, then two links in the direction of interest, and 
then back again transversely. We use a similar method for blocking by 3 and then iterate as necessary to 
get the various lattice sizes used. (We did not simulate the evolution of 96x96x48 for Nq = 15 because 
of memory limitations.) 
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FIG. 15: Linear growth of the total magnetic energy with time for several different lattice spacings, 
at fixed physical volume L x x L y x L z = (64/A: max ) x (64//c max ) X (32/A; max ), for the hard particle 
distributions (a) Afo = 7 and (b) = 15. 



N„=7; 1 , m =24,6; squeeze=3.5; fixed physical volume 

il max max * r J 

30 




FIG. 16: The curves of Fig. 15a, all shown here as dotted lines, superposed with solid lines 
corresponding to the fits of Eq. (3.9). 



for 7*t > 10. The parameters of the fit are s, t , and t\. We take the slope a (the slope of the 
tanh at zero argument) to be our result for d£^ ot /dt. The curving of the tanh is controlled 
by t\, and t± should go to infinity as we approach the continuum limit. 

Fig. 16 shows the tanh fits for the simulations of Fig. 15a: the fits work extremely well. 
The solid circles in Fig. 17 show how the fit of the slope o depends on the lattice spacing. 
The x axis is chosen to be the square of the lattice spacing because the discretization errors 
in our lattice implementation first arise at this order. Extrapolating by eye to the continuum 
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FIG. 17: Fit parameters for the (a) Nq = 7 and (b) Nq = 15 simulations of Fig. 15 as a function 
of lattice spacing squared. Solid circles are the slope a (and so d£^ ot /dt) in units of m^^/g 2 , and 
open circles are 1/ti in units of (100/7*) -1 . 



limit, we estimate that our default lattice spacing of ak m£LX = 1 for these distributions has 
lattice spacing errors no larger than roughly 10%. In contrast, the open circles in Figs. 17 
show the behavior of 1/ti, which is a lattice artifact and approaches zero in the continuum 
limit (corresponding to purely linear growth). 



IV. CONCLUSION 

The goal of this paper has been to understand how, in the weak coupling limit, the late- 
time behavior of Weibel instabilities scales with hard particle anisotropy. We can use the 
smallness of 9 = v z , characterizing the angular distribution of hard particles, as a measure 
of anisotropy. Through simulations, we have examined the slope d£^ ot /dt of the late-time 
linear growth in the total magnetic energy of soft gauge fields and found that the scaling 
of this slope is consistent with 6~ 2 and not consistent with 6~ l or 0~ 3 . If we accept the 
simple model outlined in section IB of the physics behind dS^ t /dt, this result implies that 
the limiting magnetic field strength of Weibel unstable modes scales with anisotropy as 
6~ l and is of order 

. (4.1) 
6g K } 

Of course, it would be better not to rely on such indirect arguments. A goal for future work 
should be to check the consistency of this conclusion with alternative measurements. 
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APPENDIX A: THE NIELSEN-OLESEN LIMIT 

At one time, it was conjectured that exponential instability growth would continue be- 
yond the point where non-abelian interactions became important because the fields would 
dynamically align themselves into a commuting set of color directions [22]. The fields would 
then be effectively abelian and could continue growing, just like the purely abelian case 
shown by the dashed line in Fig. 2. This conjecture seemed borne out by early simulations 
in one spatial dimension [6], such as shown by the dotted line in that figure. At the time, 
Berndt Miiller [23] predicted that three-dimensional instability growth would eventually 
have to stop because, even if the gauge fields did abelianize, Nielsen-Olesen instabilities [24] 
would eventually destroy nearly-abelian configurations as the fields continue to grow. 

In this appendix, we will discuss the largest field strength allowed for nearly-abelian 
Weibel unstable modes before Nielsen-Olesen instabilities appear. A brief review of Nielsen- 
Olesen instabilities in the context of Weibel instabilities for moderate anisotropy can be 
found in Ref. [25]. Here, we generalize to the case of extreme anisotropy. 

First imagine a situation where there was a large, constant, homogeneous magnetic field 
B that lies within an abelian subgroup of the non-abelian gauge group. In our application, 
this could represent a large magnetic field that was created by the Weibel instability and 
that might have abelianized due to non-linear dynamics, and that we are looking at this 
field on small enough time and distance scales that we can treat it as constant. For the sake 
of definiteness, consider SU(2) gauge theory and a background magnetic field 

B?(x,t) = B 5 a3 , (Al) 

where a is the adjoint color index. Now one can investigate the dispersion relation of fluctu- 
ations about this background field, including fluctuations involving other, non-commuting 
color directions. Ignoring hard particle effects, the result is [24] 

u 2 = q\ + {n + \)2\Q\gB Q - 2m s QgB , (A2) 

where q\\ is the component of momentum parallel to B ; n is the Landau orbit quantum 
number for circular motion in the plane transverse to B Q ; m s = or ±1 is the component, in 
the direction of B , of the spin of a gauge excitation; and Q = or ±1 is the charge of that 
excitation under the color generator T 3 . The last term in (A2) represents the interaction 
energy of a gauge particle's magnetic moment with the magnetic field. If we look at the 
lowest Landau orbits (n = 0) and the sector Q = m s = ±1, we get 

u 2 = qf - gB . (A3) 

uj 2 is then negative for q\\ < \/gB . This is the Nielsen-Olesen instability. 

In our application, the magnetic fields are not homogeneous. In order to make that 
approximation, the radius associated with the lowest Landau orbit should be small enough 
to fit in a region of roughly constant magnetic field. This radius is R ~ 1/y/gBo- The 
dominant Weibel-unstable modes have B roughly orthogonal to the z axis, and wavenumber 
of order (k±, k z ) ~ (m^, k max ). So, to produce the Nielsen-Olesen instability, we need both 
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R <C l/k± and R <^ l/k z . In the case of extreme anisotropy, the latter is the stricter 
constraint, equivalent to 

-= < t , (A4) 

which requires 
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Bo»-^- (A5) 
As the magnetic field grows, Nielsen-Olesen instabilities will then first appear for 

Bo - 9| • (A6) 

This corresponds to the case v = 2 in (1.17). 

There are a few approximations left to check. First, the Nielsen-Olesen analysis assumed 
that Bq was constant in time. The typical Nielsen-Olesen instability growth times generated 
by (A3) will be of order 1 / yfgB. This will be small compared to the (abelian- 

ized) Weibel instability growth time 7" 1 ~ when (A5) is satisfied. Secondly, we have 
ignored hard particle effects throughout. In general, the magnitude of the soft self-energy 
IT due to hard particles can be as large as order kl iax ~ m^/O 2 , depending on direction and 
oo/k. But this is small enough that II will be a small correction to gBo in the dispersion 
relation (A3) when (A5) is satisfied. 

The results of this paper suggest that v — 1 rather than v — 2. It may well be that 
Nielsen-Olesen effects limit the Weibel instability growth of nearly abelian fields. The results 
of this paper simply suggest that, if one starts with large amplitude, nonperturbative, non- 
abelian initial conditions, then generic non-perturbative interactions are sufficient to stop 
growth earlier, at lower field strength than (A6). 



APPENDIX B: THE COEFFICIENTS ft e 

In principle, anyone who wanted to know the specific values of fli for our distributions 
could reproduce them from the procedure outlined in the text. However, we found that 
avoiding numerical round-off errors in determining fli for large Nq required some care. So 
we will explicitly give our distributions here. 

N n = 1: Q = 1, n 2 = -1/VE 

N n = 2: Q = 1, Q 2 = -2^5/7, Q 4 = 1/7 

N n = 3: e*i = 0.585310; fi = 1, fi 2 = -0.845154, fi 4 = 0.474960, fi 6 = -0.148398 
N n = 4: «i = 0.482673; Q = 1, ft 2 = -0.907458, fi 4 = 0.589250, Q 6 = -0.277350, 
Q 8 = 0.063395 

N n = 5: c*i = 0.350109, a 2 = 0.759931; tt = 1, Q 2 = -0.971104, fi 4 = 0.730670, 
Q 6 = -0.449643, Q 8 = 0.216519, Sl w = -0.063736 

N Q = 6: ai = 0.292253, a 2 = 0.672147; fi = 1, H 2 = -0.998631, fi 4 = 0.794155, 
fi 6 = -0.544163, Q 8 = 0.311154, Q w = -0.138336, Q 12 = 0.032712 

Nq = 7: ai = 0.228412, a 2 = 0.545787, a 3 = 0.845543; Q = 1, Q 2 = -1.026410, 
Q 4 = 0.864186, Q 6 = -0.650680, Q 8 = 0.436901, Q w = -0.249727, Q 12 = 0.113554, 
Q l4 = -0.032709 
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Nq = 8: a x = 0.194863, a 2 = 0.478948, a 3 = 0.776881; fi = 1, ft 2 = -1.040916, 
tt 4 = 0.901356, Q 6 = -0.712570, Q 8 = 0.512936, Q w = -0.329365, Q 12 = 0.179614, 
n u = -0.077926, Q w = 0.018926 

N n = 9: ax = 0.159625, a 2 = 0.400532, a 3 = 0.671219, a 4 = 0.892835; tt = 1, 
Q 2 = -1.055491, Q A = 0.940673, Q 6 = -0.778554, Q 8 = 0.600618, Q 10 = -0.424975, 
n 12 = 0.271564, n u = -0.149879, Q 16 = 0.066253, Q 18 = -0.018902, 

Nq = 10: «! = 0.138840, a 2 = 0.353981, a 3 = 0.608348, a A = 0.839259; Q = 1, 
fi 2 = -1.064055, fi 4 = 0.963974, fi 6 = -0.819732, Q 8 = 0.656358, Q w = -0.490737, 
n 12 = 0.337910, tt u = -0.210147, Sl 16 = 0.111992, Q 18 = -0.047971, Q 20 = 0.011893 

Nq = 11: ai = 0.117460, a 2 = 0.303204, a 3 = 0.531249, a A = 0.753079, a 5 = 
0.921475; Q = 1, fi 2 = -1.072643, Q 4 = 0.988135, tt 6 = -0.862564, tt 8 = 0.717049, 
n w = -0.563740, n 12 = 0.417060, Sl u = -0.285605, Q w = 0.177636, Q 18 = -0.096139, 
Q 20 = 0.041819, Q 22 = -0.0118 75 

N n = 12: a x = 0.103801, a 2 = 0.270617, a 3 = 0.481526, a A = 0.697580, a 5 = 
0.879018; Q = 1, Q 2 = -1.078115, Q 4 = 1.003612, tt 6 = -0.890930, tt 8 = 0.757648, 
n w = -0.614860, Q 12 = 0.473819, Sl u = -0.343760, tt w = 0.230949, Q 18 = -0.141010, 
Q 20 = 0.074196, Q 22 = -0.031545, Q 24 = 0.007947 

Nq = 13: cti = 0.089898, a 2 = 0.236254, a 3 = 0.425611, a 4 = 0.627277, a 5 = 0.808562, 
a 6 = 0.940062; Q = 1, tt 2 = -1.083597, tt 4 = 1.019484, Q 6 = -0.920055, Q 8 = 0.800595, 
n 10 = -0.669524, Q l2 = 0.537151, Q u = -0.410252, tt w = 0.296157, Q 18 = -0.198861, 
Q 20 = 0.121655, Q 22 = -0.065052, Q 24 = 0.028013, Q 26 = -0.007936 

N n = 14: cti = 0.080480, a 2 = 0.212901, a 3 = 0.387457, a 4 = 0.579149, a 5 = 0.760545, 
a 6 = 0.905790; Q = 1, Q 2 = -1.087305, tt 4 = 1.030254, Q 6 = -0.940283, Q 8 = 0.830602, 
Q 10 = -0.708871, Q 12 = 0.583371, Sl u = -0.460883, tt w = 0.347161, Q 18 = -0.247263, 
Q 20 = 0.163617, Q 22 = -0.098747, Q 24 = 0.051563, Q 26 = -0.021818, Q 28 = 0.005568 

Nq = 15: a 1 = 0.070949, a 2 = 0.188716, a 3 = 0.346338, a 4 = 0.523712, a 5 = 0.698217, 
a 6 = 0.847595, a 7 = 0.952782; Q = 1, ft 2 = -1.091017, Q A = 1.041226, Q 6 = -0.960895, 
Q 8 = 0.861825, n 10 = -0.750085, Q 12 = 0.633157, Q 14 = -0.516209, Sl 16 = 0.405183, 
Q 18 = -0.303840, Q 20 = 0.215905, Q 22 = -0.143127, Q 2A = 0.086609, Q 26 = -0.045945, 
Q 28 = 0.019653, fi 30 = -0.005561 

Nq = 17: a x = 0.057385, a 2 = 0.153941, a 3 = 0.286217, a 4 = 0.440599, a 5 = 0.601196, 
a & = 0.751477, a 7 = 0.875973, a 8 = 0.961858; Q = 1, Q 2 = -1.096274, Q 4 = 1.056905, 
Q 6 = -0.990855, Q 8 = 0.907738, Q w = -0.812101, Q 12 = 0.709497, Q 14 = -0.603863, 
Q w = 0.499893, Q 18 = -0.400800, Q 20 = 0.310023, Q 22 = -0.229608, Q 24 = 0.161402, 
Q 26 = -0.106052, Q 28 = 0.063689, Q 30 = -0.033601, Q 32 = 0.014305, Q 34 = -0.004046 

Nq = 19: ai = 0.047353, a 2 = 0.127821, a 3 = 0.239919, a 4 = 0.374139, a 5 = 0.519102, 
a 6 = 0.662513, a 7 = 0.792211, a 8 = 0.897196, a 9 = 0.968556; tt = 1, Q 2 = -1.100133, 
Q 4 = 1.068574, Q 6 = -1.013442, Q 8 = 0.942914, Q w = -0.860525, Q l2 = 0.770493, 
n 1A = -0.675846, Q w = 0.580283, Q 18 = -0.486451, Q 20 = 0.397320, Q 22 = -0.314853, 
VL 2A = 0.241036, Q 26 = -0.176940, Q 28 = 0.123418, Q 30 = -0.080579, Q 32 = 0.048126, 
Q u = -0.025290, Q 36 = 0.010730, Q 38 = -0.003034 

Nq = 21: oil = 0.039729, a 2 = 0.107746, a 3 = 0.203686, a 4 = 0.320731, a 5 = 0.450563, 
a & = 0.583957, a 7 = 0.711432, a 8 = 0.823930, a 9 = 0.913455, a w = 0.973637; Q = 1, 
Q 2 = -1.103050, fi 4 = 1.077489, Q 6 = -1.030870, Q 8 = 0.970394, Q 10 = -0.898899, 
Q 12 = 0.819660, Q 14 = -0.735041, tt w = 0.647973, Q 18 = -0.560615, Q 20 = 0.475461, 
Q 22 = -0.394267, Q 2A = 0.318906, Q 26 = -0.250564, Q 28 = 0.190366, Q 30 = -0.138832, 
Q 32 = 0.096280, Q 34 = -0.062563, Q 36 = 0.037213, Q 38 = -0.019497, Q 40 = 0.008251, 
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n 42 = -0.002333 

Nq = 23: ai = 0.033803, a 2 = 0.092008, a 3 = 0.174892, a 4 = 0.277451, a 5 = 0.393487, 
a 6 = 0.515992, a 7 = 0.637566, a 8 = 0.750864, a 9 = 0.849043, a 10 = 0.926171, a u = 
0.977582; Q = 1, Q 2 = -1.105308, Q 4 = 1.084451, Q 6 = -1.044589, Q 8 = 0.992235, 
Qio = -0.929737, Q 12 = 0.859695, Sl u = -0.783977, Q w = 0.704930, Q 18 = -0.624313, 
Q 20 = 0.544203, Q 22 = -0.466119, Vt 2A = 0.391741, Q 26 = -0.322219, Q 28 = 0.258733, 
Q 30 = -0.201979, Q 32 = 0.152571, Q u = -0.110714, Q 36 = 0.076443, Q 38 = -0.049493, 
fi 40 = 0.029346, tt A2 = -0.015341, fi 44 = 0.006479, fi 46 = -0.001832 

Nq = 25: «i = 0.029107, a 2 = 0.079452, a 3 = 0.151682, a 4 = 0.242041, a 5 = 0.345834, 
a 6 = 0.457666, a 7 = 0.571724, a 8 = 0.682079, a 9 = 0.782997, a w = 0.869230, a u = 
0.936297, a 12 = 0.980705; Q = 1, Sl 2 = -1.107092, Q 4 = 1.089990, tt 6 = -1.055575, 
Q 8 = 1.009862, Q w = -0.954847, Q 12 = 0.892630, Q 14 = -0.824713, Q w = 0.752995, 
Q 18 = -0.678916, Q 20 = 0.604201, Q 22 = -0.530149, Vt 2A = 0.458227, Q 26 = -0.389494, 
Q 28 = 0.325078, Q 30 = -0.265727, Q 32 = 0.212172, Q 34 = -0.164805, Q 36 = 0.123933, 
Q 38 = -0.089584, Q 40 = 0.061641, Q 42 = -0.039796, Q 44 = 0.023538, Q m = -0.012284, 
Q 48 = 0.005180, Q 50 = -0.001465 

Various perturbative results for instabilities can be calculated directly from the f2/'s using 
the following formula for the transverse gluon self-energy in the special case that the gluon 
momentum k points along the beam axis [12]: 



U±(u, ke z 
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with 



K e (rj) = (1 + rf)5 m + (1 - V 2 )[(£ + l)Q e+1 ( V ) - (I - l)vQi(v)}- 



(Bib) 



Here, Qi(t]) is the Legendre function of the second kind defined so that it is regular at rj = oo 
and the cut is chosen to run from —1 to +1. For example, Qo{z) = |ln[(z + l)/(z — 1)]. 
The corresponding dispersion relation is 



— id 2 + k 2 + IT_i_(u;, ke 2 



0. 



(B2) 



For a given distribution Q(9), one can solve this equation numerically for each k, which 
is how Fig. 7 was produced. By scanning over k, the largest growth rate 7 = Imu and 
corresponding wavenumber k* can be found. The remaining parameters in Table II are 
given in terms of the Qi as 



and [2] 
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APPENDIX C: DAMPING LARGE I AND m MODES 



In previous work for moderate anisotropy [13], we found we could reduce errors from a 
finite £ max cut-off by damping t modes near the cut-off. There, we modified the equations 
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of motion for the W^ m 's to 



= (original) - 7dam P W^ m Q{\ + £- £ damp ) , (CI) 

where 6(z) is the step function. This introduced damping for all modes with £ between 
4am P and £ max . We chose 

^damp |_3^maxJ i 7damp 1 ,. (^^) 

V *-max 

The rationale behind the choice of 7dam P is discussed in Ref. [13]. In brief, energy in W^ m 's 
of large £, m tends to cascade to still higher £, m. The presence of a cutoff can "reflect" the 
energy back to low £, m, which is unphysical. Damping avoids this by absorbing this energy, 
which reproduces the physics of its cascading to arbitrarily high £, m. 

In the current work, we introduce similar damping near the cut-off m max on m. Also, we 
have changed the procedure to turn on the amount of damping more gradually as £ or m 
increase. 13 This prevents reflection at the boundary between modes which are and are not 
damped. In this paper, we replace (CI) by 



dWjm 

dt 



= (original) - 7* m W/ m 6(| + £ - £ damp ) 6(± + m - m damp ), (C3) 



with 



^ _ Moo (|+l-ldamp) | (| + 777 - 777 damp ) 

\Z^max (^max ^damp) ^/77l max (77l max 77l damp ) 

^damp |_3^maxJ and 77l damp |_377l max J . (C5) 
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